=> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)
J’indique un rayon en mètres (< 15000) pour sélectionner les parcelles autour du point
Code
# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")coord_gmaps<-"43.447894436406216, 1.2886163291688764"lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify =TRUE)[,1])lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify =TRUE)[,2])# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))rayon<-10000# création d'une table sf «point» avec les coordonnées saisies# transformation des coordonnées en syst de proj 2154 (Lambert II - Français) point<-data.frame(lon,lat,rayon) %>%st_as_sf(coords =c("lon","lat"),crs ="EPSG:4326") %>%mutate(coord_pt_gps=st_as_text(geometry)) %>%st_transform("EPSG:2154") %>%st_sf() %>%clean_names() %>%rename(geom=geometry)# st_crs(point)# st_geometry_type(point)# plot(point)
La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.
1.2 Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me
Code
# le mot de passe est stocké dans un secret Vault# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")# postgresql_password <- "1tfawt3nj7fgzo3w7cma"# Sys.setenv("PASS_POSTGRESQL"="1tfawt3nj7fgzo3w7cma")# Connection à PostgreSQL cnx <-dbConnect(PostgreSQL(),user ="projet-funathon",password =Sys.getenv("PASS_POSTGRESQL"),host ="postgresql-758156",dbname ="defaultdb",port =5432,options="-c search_path=rpg,public") # specify what schema to connect to
1.3 je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.
Code
# suppression de la table «point» si elle existedbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")
<PostgreSQLResult>
Code
# écriture de la table point dans la base PostGiswrite_sf(point, cnx, append = T)# ajout d'une clé primairedbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")
<PostgreSQLResult>
Code
# ajout d'un indexdbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")
<PostgreSQLResult>
Code
# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis -------dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")
<PostgreSQLResult>
Code
dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.* FROM rpg.point p, rpg.parcelles r WHERE ST_DWithin(p.geom,r.geom,p.rayon);")
# ajout d'un indexdbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")
<PostgreSQLResult>
Code
# 4 - Téléchargement des parcelles proches sous R-------------------------------# parc_prox<-read_sf(cnx,parc_prox)parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")
1.4 J’affiche les parcelles autour de mon point avec une carte interactive leaflet
Code
#plot(st_geometry(parc_prox))#plot(st_geometry(point),add = T,col = "red")# Transformation de la projection car leaflet ne connait que le WGS 84carte_parc_prox<- parc_prox %>%st_transform(4326)# Marqueur du pointpt_mark<- point %>%st_transform(4326)# ajout du libellé des culturescarte_parc_prox_lib<-carte_parc_prox %>%left_join(lib_cult %>%select(-code_groupe_culture),by=c("code_cultu"="code_culture")) #création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.labels<-sprintf("<strong>id_parcel : </strong>%s<br/> <strong>Groupe culture : </strong>%s<br/> <strong>Culture : </strong>%s<br/> <strong>Surface (ha) : </strong>%s<br/> <strong>Département : </strong>%s<br/> <strong>Commune : </strong>%s<br/>", parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture, carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc, parc_prox$insee_dep,parc_prox$nom_com) %>%lapply(htmltools::HTML)# labels# création d'une palette de couleurs associée au groupe de culturefactpal <-colorFactor("Paired", parc_prox$code_group)carte_parc_prox_html <-leaflet(carte_parc_prox_lib) %>%# addProviderTiles("Esri.WorldImagery") %>%# addTiles() %>% addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%addPolygons( #fillColor="white",fillColor=~factpal(code_group),weight =2,opacity =1,color ="#ffd966",dashArray ="3",fillOpacity =0.5,highlight =highlightOptions(weight =5,color ="#A40000",dashArray ="",fillOpacity =0.0,bringToFront =TRUE),label = labels,labelOptions =labelOptions(style =list("font-weight"="normal", padding ="3px 8px"),textsize ="15px",direction ="auto",encoding="UTF-8")) %>%addMarkers(data=pt_mark,~lon,~lat, popup =~coord_pt_gps, label =~coord_pt_gps)
Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Paired is 12
Returning the palette you asked for with that many colors
Storing counts in `nn`, as `n` already present in input
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
t1 %>%setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>%kable(format="html",caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(t1), bold = T, color ="white", background ="grey")
Structure des cultures au niveau local
code
groupe de cultures
nombre de parcelles
(%)
surface (ha)
surface (%)
taille moyenne (ha)
1
Blé tendre
631
10,9
3 044,52
17,6
4,8
2
Maïs grain et ensilage
368
6,4
2 922,39
16,9
7,9
4
Autres céréales
404
7,0
1 826,24
10,6
4,5
11
Gel (surfaces gelées sans production)
1 437
24,9
1 494,91
8,7
1,0
6
Tournesol
308
5,3
1 430,40
8,3
4,6
18
Prairies permanentes
573
9,9
1 260,24
7,3
2,2
7
Autres oléagineux
184
3,2
1 237,37
7,2
6,7
3
Orge
214
3,7
947,61
5,5
4,4
16
Fourrage
234
4,1
804,31
4,7
3,4
19
Prairies temporaires
302
5,2
802,74
4,6
2,7
8
Protéagineux
116
2,0
538,13
3,1
4,6
5
Colza
106
1,8
468,59
2,7
4,4
28
Divers
754
13,1
184,22
1,1
0,2
15
Légumineuses à grains
34
0,6
147,28
0,9
4,3
25
Légumes ou fleurs
38
0,7
46,95
0,3
1,2
21
Vignes
21
0,4
37,31
0,2
1,8
17
Estives et landes
16
0,3
33,19
0,2
2,1
24
Autres cultures industrielles
23
0,4
32,70
0,2
1,4
20
Vergers
7
0,1
4,76
0,0
0,7
Total
-
5 770
100,0
17 263,86
100,1
3,0
Code
# rm(t1)
1.6 Comparaison avec la répartition des cultures au niveau départemental et national
Code
# jointure spatiale avec la couche département pour récupérer le département où tombe le point# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)dep <-s3read_using(FUN = sf::read_sf,layer ="departement",object ="2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",bucket ="projet-funathon",opts =list("region"="")) %>%st_transform(2154)# jointuredf<-point %>%st_join(dep) %>%st_drop_geometry() %>%select(insee_dep)dep_pt<-df[1,1]rm(df)# sinon sélection du département regroupant la plus grande surface agricole# cas où le cercle recouvre plusieurs départements# df<-parc_prox %>% st_drop_geometry() %>% # count(insee_dep,wt=surf_parc) %>% # arrange(desc(n))# dep_pt<-df[1,1]# rm(df)# calcul des % surfaces autour du pointstat_pt <- parc_prox %>%st_drop_geometry() %>%count(code_group,wt=surf_parc) %>%add_tally(n) %>%mutate(pct_surf_local=round(100*n/nn,1)) %>%select(code_group, pct_surf_local)
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
# récup des % surfaces départementalesstat_dep_pt<-s3read_using(FUN=readr::read_rds,object ="2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",bucket ="projet-funathon",opts =list("region"="")) %>%filter(insee_dep %in% dep_pt) %>%select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>%rename(pct_surf_dep = pct_surf)# récup des % surfaces nationalesstat_fm<-s3read_using(FUN=readr::read_csv,object ="2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",col_types =cols(code_group =col_character()),bucket ="projet-funathon",opts =list("region"="")) %>%select(code_group,libelle_groupe_culture,pct_surf) %>%rename(pct_surf_fm = pct_surf)# appariement des stas locales, départementales, nationalesstat_compar<-stat_fm %>%left_join(stat_dep_pt %>%select(code_group, pct_surf_dep),by="code_group") %>%left_join(stat_pt,by="code_group") %>%select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>%arrange(desc(pct_surf_local)) %>%adorn_totals() stat_compar %>%setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>%kable(format="html",caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(stat_compar), bold = T, color ="white", background ="grey")
Structure des cultures
Groupe de cultures
surf. locales (%)
surf. départ. (%)
surf. France m. (%)
Blé tendre
17,6
14,3
17,7
Maïs grain et ensilage
16,9
8,6
10,0
Autres céréales
10,6
12,4
3,9
Gel (surfaces gelées sans production)
8,7
4,4
1,6
Tournesol
8,3
12,6
2,5
Prairies permanentes
7,3
17,9
27,7
Autres oléagineux
7,2
3,1
0,7
Orge
5,5
3,5
6,2
Fourrage
4,7
5,9
3,7
Prairies temporaires
4,6
4,2
5,3
Protéagineux
3,1
1,8
1,2
Colza
2,7
1,8
3,5
Divers
1,1
1,3
1,3
Légumineuses à grains
0,9
0,3
0,2
Légumes ou fleurs
0,3
0,3
1,5
Estives et landes
0,2
6,9
8,1
Vignes
0,2
0,4
2,2
Autres cultures industrielles
0,2
0,1
1,8
Vergers
0,0
0,1
0,4
Plantes à fibres
NA
0,0
0,4
Riz
NA
NA
0,0
Fruits à coque
NA
0,0
0,2
Oliviers
NA
0,0
0,0
Total
100,1
99,9
100,1
1.7 Graphique de comparaison des cultures au niveau local, départemental et national
Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local tab<-stat_compar %>%filter(libelle_groupe_culture!="Total") %>%slice_head(n=10) %>%rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)tab_piv<-tab %>%pivot_longer(!libelle_groupe_culture) %>%rename(secteur=name) # je mets à 0 les valeurs manquantes#tab_piv<-tab_piv %>% coalesce(value,0L) tab_piv[is.na(tab_piv)] <-0# levels(as.factor(tab_piv$secteur))# je réordonne les secteurs dans le "bon" ordre, avec factor tab_piv$secteur<-factor(tab_piv$secteur,levels=c("france","departement","local"))tab_piv<-tab_piv %>%arrange(desc(secteur), desc(value))# je réordonne les cultures par surface décroissante au niveau local, avec factorx <- tab_piv %>%filter(secteur =="local") %>%arrange(value) %>%select(libelle_groupe_culture)y <-pull(x, libelle_groupe_culture)tab_piv$libelle_groupe_culture <-factor(tab_piv$libelle_groupe_culture, levels = y)# ggplot => problème de tri : affichage des prairies avant le maïs ? p<-ggplot(tab_piv, aes(x =libelle_groupe_culture,y = value, fill=factor(secteur, levels=c("france","departement","local")))) +geom_col(position ="dodge") +labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y ="%", fill ="Secteur") +theme_classic()# je flippe le graphique pour avoir des barres horizontales p+coord_flip()
différences de structure locale, departementale, nationale
1.8 Je teste un graphique par secteur avec facet_wrap
Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local tab<-stat_compar %>%filter(libelle_groupe_culture!="Total") %>%slice_head(n=10) %>%rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)tab_piv<-tab %>%pivot_longer(!libelle_groupe_culture) %>%rename(secteur=name) # je mets à 0 les valeurs manquantes#tab_piv<-tab_piv %>% coalesce(value,0L) tab_piv[is.na(tab_piv)] <-0# ggplot => problème de tri : affichage des prairies avant le maïs ? ggplot(tab_piv, aes(x =reorder(libelle_groupe_culture,+value), y = value)) +geom_col(fill ="lightblue", colour ="black",position ="dodge") +labs(title="Surface par culture", x="Culture", y ="%", fill ="Secteur") +geom_text(aes(label = value), hjust =-0.3, size =8/.pt, colour ="black") +theme_classic() +coord_flip() +facet_wrap(~secteur,nrow=3,scales='free')
différences de structure locale, departementale, nationale
1.9 Pour aller plus loin, comparer la taille moyenne des parcelles selon les secteurs
Source Code
---title: "Sujet2-Etape1"format: htmlcode-fold: truecode-tools: truecode-summary: "voir le code"# embed-resources: true# standalone: trueself-contained: truetoc: truetoc-depth: 3number-sections: truetoc-location: leftcss: styles.csstheme: unitededitor: visualoutput-file: Sujet2-Etape1.html---```{r sys.setenv, include=FALSE}#| echo: false#install.packages("aws.s3", repos = "https://cloud.R-project.org")Sys.setenv("AWS_ACCESS_KEY_ID"="Y7E0FP1FZ22FQ3YP44UM","AWS_SECRET_ACCESS_KEY"="J4+4P5HnlFbvTCgR+dIhyYQWfeymzaAhM2KRyVhj","AWS_DEFAULT_REGION"="us-east-1","AWS_SESSION_TOKEN"="eyJhbGciOiJIUzUxMiIsInR5cCI6IkpXVCJ9.eyJhY2Nlc3NLZXkiOiJZN0UwRlAxRloyMkZRM1lQNDRVTSIsImFsbG93ZWQtb3JpZ2lucyI6WyIqIl0sImF1ZCI6WyJtaW5pby1kYXRhbm9kZSIsIm9ueXhpYSIsImFjY291bnQiXSwiYXV0aF90aW1lIjoxNjg2ODM3OTgxLCJhenAiOiJvbnl4aWEiLCJlbWFpbCI6ImJlcnRyYW5kLmJhbGxldEBhZ3JpY3VsdHVyZS5nb3V2LmZyIiwiZW1haWxfdmVyaWZpZWQiOnRydWUsImV4cCI6MTY4NjkyNDM5NiwiZmFtaWx5X25hbWUiOiJCQUxMRVQiLCJnaXZlbl9uYW1lIjoiQmVydHJhbmQiLCJncm91cHMiOlsiZnVuYXRob24iXSwiaWF0IjoxNjg2ODM3OTgyLCJpc3MiOiJodHRwczovL2F1dGgubGFiLnNzcGNsb3VkLmZyL2F1dGgvcmVhbG1zL3NzcGNsb3VkIiwianRpIjoiNjJhYjQ5ZWEtYmU0Yi00MzE4LTljYjktMWY1NGZkNmRmYTkwIiwibmFtZSI6IkJlcnRyYW5kIEJBTExFVCIsIm5vbmNlIjoiODllYmM1YzYtNWZiYy00ODRkLTk3MDEtOTJkMTY0NzJjNjFjIiwicG9saWN5Ijoic3Rzb25seSIsInByZWZlcnJlZF91c2VybmFtZSI6ImJiYWxsZXQiLCJyZWFsbV9hY2Nlc3MiOnsicm9sZXMiOlsib2ZmbGluZV9hY2Nlc3MiLCJ1bWFfYXV0aG9yaXphdGlvbiJdfSwicmVzb3VyY2VfYWNjZXNzIjp7ImFjY291bnQiOnsicm9sZXMiOlsibWFuYWdlLWFjY291bnQiLCJtYW5hZ2UtYWNjb3VudC1saW5rcyIsInZpZXctcHJvZmlsZSJdfX0sInNjb3BlIjoib3BlbmlkIHByb2ZpbGUgZ3JvdXBzIGVtYWlsIiwic2Vzc2lvblBvbGljeSI6ImV5SldaWEp6YVc5dUlqb2lNakF4TWkweE1DMHhOeUlzSWxOMFlYUmxiV1Z1ZENJNlczc2lSV1ptWldOMElqb2lRV3hzYjNjaUxDSkJZM1JwYjI0aU9sc2ljek02S2lKZExDSlNaWE52ZFhKalpTSTZXeUpoY200NllYZHpPbk16T2pvNmNISnZhbVYwTFdaMWJtRjBhRzl1SWl3aVlYSnVPbUYzY3pwek16bzZPbkJ5YjJwbGRDMW1kVzVoZEdodmJpOHFJbDE5TEhzaVJXWm1aV04wSWpvaVFXeHNiM2NpTENKQlkzUnBiMjRpT2xzaWN6TTZUR2x6ZEVKMVkydGxkQ0pkTENKU1pYTnZkWEpqWlNJNld5SmhjbTQ2WVhkek9uTXpPam82S2lKZExDSkRiMjVrYVhScGIyNGlPbnNpVTNSeWFXNW5UR2xyWlNJNmV5SnpNenB3Y21WbWFYZ2lPaUprYVdabWRYTnBiMjR2S2lKOWZYMHNleUpGWm1abFkzUWlPaUpCYkd4dmR5SXNJa0ZqZEdsdmJpSTZXeUp6TXpwSFpYUlBZbXBsWTNRaVhTd2lVbVZ6YjNWeVkyVWlPbHNpWVhKdU9tRjNjenB6TXpvNk9pb3ZaR2xtWm5WemFXOXVMeW9pWFgxZGZRPT0iLCJzZXNzaW9uX3N0YXRlIjoiMzAwYzgzNDktZWQ5MS00MGVhLThhZGQtYTkxYzlhMTE5YzlkIiwic2lkIjoiMzAwYzgzNDktZWQ5MS00MGVhLThhZGQtYTkxYzlhMTE5YzlkIiwic3ViIjoiNGQ1N2NlN2QtOWNhYS00MDk4LWFjNDktNzY0MjAyYzJhN2I5IiwidHlwIjoiQmVhcmVyIn0.rrSRJ9eWjRyk-SPG-1-UIyLK4EsS40ViBwstKwrSim20xEZQb9f9-jEPNBeajvGc4k3N-zGLQeyHWu7IvfkKKQ","AWS_S3_ENDPOINT"="minio.lab.sspcloud.fr")``````{r setup, include=FALSE}#| label: setup#| echo: falselibrary(tidyverse) library(aws.s3)library(sf)library(RPostgreSQL)library(janitor)if (!"kableExtra"%in%installed.packages()) { install.packages("kableExtra") }library(kableExtra)library(leaflet)library(htmlwidgets)library(RColorBrewer)# install.packages("aws.s3", repos = "https://cloud.R-project.org")bucketlist(region="")``````{r acces donnees, include=FALSE}#| echo: false# données -----------------------------------------------------------------aws.s3::get_bucket("projet-funathon", region ="", prefix ="2023/sujet2/diffusion/ign/rpg")``````{r, include=FALSE }#| echo: false#| label: lecture libellés groupes cultureslib_cult<-s3read_using(FUN = read_csv2, object ="2023/sujet2/diffusion/ign/rpg/REF_CULTURES_GROUPES_CULTURES_2020.csv",col_types =cols(.default =col_character()),bucket ="projet-funathon",opts =list("region"="")) %>%clean_names() lib_group_cult<-lib_cult %>%select(code_groupe_culture,libelle_groupe_culture) %>%distinct(code_groupe_culture,.keep_all=T)```## Etape 1 : Première manipulation du RPG### Je choisis un point sur la carte : <https://www.google.fr/maps>=\> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)J'indique un rayon en mètres (\< 15000) pour sélectionner les parcelles autour du point```{r}#| label: choix coordonnées + rayon# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")coord_gmaps<-"43.447894436406216, 1.2886163291688764"lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify =TRUE)[,1])lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify =TRUE)[,2])# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))rayon<-10000# création d'une table sf «point» avec les coordonnées saisies# transformation des coordonnées en syst de proj 2154 (Lambert II - Français) point<-data.frame(lon,lat,rayon) %>%st_as_sf(coords =c("lon","lat"),crs ="EPSG:4326") %>%mutate(coord_pt_gps=st_as_text(geometry)) %>%st_transform("EPSG:2154") %>%st_sf() %>%clean_names() %>%rename(geom=geometry)# st_crs(point)# st_geometry_type(point)# plot(point)```La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.### Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me```{r}#| label: connection PostGis# le mot de passe est stocké dans un secret Vault# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")# postgresql_password <- "1tfawt3nj7fgzo3w7cma"# Sys.setenv("PASS_POSTGRESQL"="1tfawt3nj7fgzo3w7cma")# Connection à PostgreSQL cnx <-dbConnect(PostgreSQL(),user ="projet-funathon",password =Sys.getenv("PASS_POSTGRESQL"),host ="postgresql-758156",dbname ="defaultdb",port =5432,options="-c search_path=rpg,public") # specify what schema to connect to```### je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.```{r}#| label: requete SQL selection parcelles# suppression de la table «point» si elle existedbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")# écriture de la table point dans la base PostGiswrite_sf(point, cnx, append = T)# ajout d'une clé primairedbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")# ajout d'un indexdbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis -------dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.* FROM rpg.point p, rpg.parcelles r WHERE ST_DWithin(p.geom,r.geom,p.rayon);")# ajout d'une clé primairedbSendQuery(cnx,"ALTER TABLE rpg.parc_prox ADD CONSTRAINT parc_prox_pk PRIMARY KEY(id_parcel);")# ajout d'un indexdbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")# 4 - Téléchargement des parcelles proches sous R-------------------------------# parc_prox<-read_sf(cnx,parc_prox)parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")```### J'affiche les parcelles autour de mon point avec une carte interactive leaflet```{r}#| label: affichage parcelles#plot(st_geometry(parc_prox))#plot(st_geometry(point),add = T,col = "red")# Transformation de la projection car leaflet ne connait que le WGS 84carte_parc_prox<- parc_prox %>%st_transform(4326)# Marqueur du pointpt_mark<- point %>%st_transform(4326)# ajout du libellé des culturescarte_parc_prox_lib<-carte_parc_prox %>%left_join(lib_cult %>%select(-code_groupe_culture),by=c("code_cultu"="code_culture")) #création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.labels<-sprintf("<strong>id_parcel : </strong>%s<br/> <strong>Groupe culture : </strong>%s<br/> <strong>Culture : </strong>%s<br/> <strong>Surface (ha) : </strong>%s<br/> <strong>Département : </strong>%s<br/> <strong>Commune : </strong>%s<br/>", parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture, carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc, parc_prox$insee_dep,parc_prox$nom_com) %>%lapply(htmltools::HTML)# labels# création d'une palette de couleurs associée au groupe de culturefactpal <-colorFactor("Paired", parc_prox$code_group)carte_parc_prox_html <-leaflet(carte_parc_prox_lib) %>%# addProviderTiles("Esri.WorldImagery") %>%# addTiles() %>% addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%addPolygons( #fillColor="white",fillColor=~factpal(code_group),weight =2,opacity =1,color ="#ffd966",dashArray ="3",fillOpacity =0.5,highlight =highlightOptions(weight =5,color ="#A40000",dashArray ="",fillOpacity =0.0,bringToFront =TRUE),label = labels,labelOptions =labelOptions(style =list("font-weight"="normal", padding ="3px 8px"),textsize ="15px",direction ="auto",encoding="UTF-8")) %>%addMarkers(data=pt_mark,~lon,~lat, popup =~coord_pt_gps, label =~coord_pt_gps)# leaflet(pt_mark) %>% addTiles() %>% addMarkers(~lon,~lat)carte_parc_prox_html# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox_html.html")```### Quelle est la structure des parcelles agricoles autour de mon point ?```{r}#| label: stats sur les groupes de cultures #| tbl-cap: Structure des cultures au niveau local t1 <- parc_prox %>%st_drop_geometry() %>%count(code_group) %>%add_tally(n) %>%mutate(n_pct=round(100*n/nn,1)) %>%select(-nn) %>%rename(n_parcelles=n) %>%# adorn_totals() %>% cbind(# comptage des surfacesparc_prox %>%st_drop_geometry() %>%count(code_group,wt=surf_parc) %>%add_tally(n) %>%mutate(surf_pct=round(100*n/nn,1)) %>%select(-nn) %>%rename(surf_parc_ha=n) %>%select(surf_parc_ha, surf_pct) # %>% # adorn_totals()) %>%left_join(lib_group_cult,by=c("code_group"="code_groupe_culture")) %>%select(code_group,libelle_groupe_culture,everything()) %>%# arrange(as.numeric(code_group)) %>% arrange(desc(surf_parc_ha)) %>%adorn_totals() %>%mutate(taille_moy_parc=round(surf_parc_ha/n_parcelles,1))t1 %>%setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>%kable(format="html",caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(t1), bold = T, color ="white", background ="grey")# rm(t1)```### Comparaison avec la répartition des cultures au niveau départemental et national```{r}#| label: comparaison structure départementale et nationale #| tbl-cap: Structure des cultures # jointure spatiale avec la couche département pour récupérer le département où tombe le point# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)dep <-s3read_using(FUN = sf::read_sf,layer ="departement",object ="2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",bucket ="projet-funathon",opts =list("region"="")) %>%st_transform(2154)# jointuredf<-point %>%st_join(dep) %>%st_drop_geometry() %>%select(insee_dep)dep_pt<-df[1,1]rm(df)# sinon sélection du département regroupant la plus grande surface agricole# cas où le cercle recouvre plusieurs départements# df<-parc_prox %>% st_drop_geometry() %>% # count(insee_dep,wt=surf_parc) %>% # arrange(desc(n))# dep_pt<-df[1,1]# rm(df)# calcul des % surfaces autour du pointstat_pt <- parc_prox %>%st_drop_geometry() %>%count(code_group,wt=surf_parc) %>%add_tally(n) %>%mutate(pct_surf_local=round(100*n/nn,1)) %>%select(code_group, pct_surf_local) # récup des % surfaces départementalesstat_dep_pt<-s3read_using(FUN=readr::read_rds,object ="2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",bucket ="projet-funathon",opts =list("region"="")) %>%filter(insee_dep %in% dep_pt) %>%select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>%rename(pct_surf_dep = pct_surf)# récup des % surfaces nationalesstat_fm<-s3read_using(FUN=readr::read_csv,object ="2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",col_types =cols(code_group =col_character()),bucket ="projet-funathon",opts =list("region"="")) %>%select(code_group,libelle_groupe_culture,pct_surf) %>%rename(pct_surf_fm = pct_surf)# appariement des stas locales, départementales, nationalesstat_compar<-stat_fm %>%left_join(stat_dep_pt %>%select(code_group, pct_surf_dep),by="code_group") %>%left_join(stat_pt,by="code_group") %>%select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>%arrange(desc(pct_surf_local)) %>%adorn_totals() stat_compar %>%setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>%kable(format="html",caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(stat_compar), bold = T, color ="white", background ="grey")```### Graphique de comparaison des cultures au niveau local, départemental et national```{r}#| label: faire un graphique comparant les structures#| fig-cap: différences de structure locale, departementale, nationale # je sélectionne les 10 groupes de cultures les plus répandus au niveau local tab<-stat_compar %>%filter(libelle_groupe_culture!="Total") %>%slice_head(n=10) %>%rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)tab_piv<-tab %>%pivot_longer(!libelle_groupe_culture) %>%rename(secteur=name) # je mets à 0 les valeurs manquantes#tab_piv<-tab_piv %>% coalesce(value,0L) tab_piv[is.na(tab_piv)] <-0# levels(as.factor(tab_piv$secteur))# je réordonne les secteurs dans le "bon" ordre, avec factor tab_piv$secteur<-factor(tab_piv$secteur,levels=c("france","departement","local"))tab_piv<-tab_piv %>%arrange(desc(secteur), desc(value))# je réordonne les cultures par surface décroissante au niveau local, avec factorx <- tab_piv %>%filter(secteur =="local") %>%arrange(value) %>%select(libelle_groupe_culture)y <-pull(x, libelle_groupe_culture)tab_piv$libelle_groupe_culture <-factor(tab_piv$libelle_groupe_culture, levels = y)# ggplot => problème de tri : affichage des prairies avant le maïs ? p<-ggplot(tab_piv, aes(x =libelle_groupe_culture,y = value, fill=factor(secteur, levels=c("france","departement","local")))) +geom_col(position ="dodge") +labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y ="%", fill ="Secteur") +theme_classic()# je flippe le graphique pour avoir des barres horizontales p+coord_flip()```### Je teste un graphique par secteur avec facet_wrap```{r}#| label: faire un graphique facet_wrap #| fig-cap: différences de structure locale, departementale, nationale # je sélectionne les 10 groupes de cultures les plus répandus au niveau local tab<-stat_compar %>%filter(libelle_groupe_culture!="Total") %>%slice_head(n=10) %>%rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)tab_piv<-tab %>%pivot_longer(!libelle_groupe_culture) %>%rename(secteur=name) # je mets à 0 les valeurs manquantes#tab_piv<-tab_piv %>% coalesce(value,0L) tab_piv[is.na(tab_piv)] <-0# ggplot => problème de tri : affichage des prairies avant le maïs ? ggplot(tab_piv, aes(x =reorder(libelle_groupe_culture,+value), y = value)) +geom_col(fill ="lightblue", colour ="black",position ="dodge") +labs(title="Surface par culture", x="Culture", y ="%", fill ="Secteur") +geom_text(aes(label = value), hjust =-0.3, size =8/.pt, colour ="black") +theme_classic() +coord_flip() +facet_wrap(~secteur,nrow=3,scales='free')```### Pour aller plus loin, comparer la taille moyenne des parcelles selon les secteurs